/*=========================================================================
 *
 *  Copyright SINAPSE: Scalable Informatics for Neuroscience, Processing and Software Engineering
 *            The University of Iowa
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
/*=========================================================================
 *
 *  Program:   Insight Segmentation & Registration Toolkit
 *  Module:    $RCSfile$
 *  Language:  C++
 *
 *  Copyright (c) Insight Software Consortium. All rights reserved.
 *  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
 *
 *  This software is distributed WITHOUT ANY WARRANTY; without even
 *  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 *  PURPOSE.  See the above copyright notices for more information.
 *
 *  =========================================================================*/
#ifndef __itkInverseConsistentLandmarks_hxx
#define __itkInverseConsistentLandmarks_hxx
#include "itkInverseConsistentLandmarks.h"
#include <cstdio>
#include <ctime>
#include <string>
#include <iostream>
#include "itkNumberToString.h"

namespace itk
{
template <typename PointStorageType, typename PointSetType>
InverseConsistentLandmarks<PointStorageType, PointSetType>::InverseConsistentLandmarks(void)
{
  ImageDims[0] = 12345;
  ImageDims[1] = 12345;
  ImageDims[2] = 12345;
  ImageDims[3] = 12345;
}

template <typename PointStorageType, typename PointSetType>
InverseConsistentLandmarks<PointStorageType, PointSetType>::InverseConsistentLandmarks(const int XDim,
                                                                                       const int YDim,
                                                                                       const int ZDim,
                                                                                       const int TDim)
{
  ImageDims[0] = XDim;
  ImageDims[1] = YDim;
  ImageDims[2] = ZDim;
  ImageDims[3] = TDim;
}

template <typename PointStorageType, typename PointSetType>
InverseConsistentLandmarks<PointStorageType, PointSetType> &
InverseConsistentLandmarks<PointStorageType, PointSetType>::operator=(const InverseConsistentLandmarks & rhs)
{
  this->std::map<std::string, PointType>::operator=(rhs);
  ImageDims[0] = rhs.ImageDims[0];
  ImageDims[1] = rhs.ImageDims[1];
  ImageDims[2] = rhs.ImageDims[2];
  ImageDims[3] = rhs.ImageDims[3];
  return *this;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadPointTypes(const std::string lmrkfilename)
{
  const Landmark_File_Format lmkff = QueryLandmarkFile(lmrkfilename);
  bool                       check = false;

  switch (lmkff)
  {
    case GEC_LANDMARKS:
    {
      check = InverseConsistentLandmarks::ReadGECPointTypes(lmrkfilename);
    }
    break;
    case INTELLX_LANDMARKS:
    {
      check = InverseConsistentLandmarks::ReadIntellXPointTypes(lmrkfilename);
    }
    break;
    case ANALYZE_LANDMARKS:
    {
      check = InverseConsistentLandmarks::ReadAnalyzePointTypes(lmrkfilename);
    }
    break;
    case IPL_LANDMARKS:
    {
      check = InverseConsistentLandmarks::ReadIPLPointTypes(lmrkfilename);
    }
    break;
    case IPL_TALAIRACH_BOUNDS_LANDMARKS:
    {
      check = InverseConsistentLandmarks::ReadIPLTalairachPointTypes(lmrkfilename);
    }
    break;
    case IPL_CEREBELLAR_BOUNDS_LANDMARKS:
    {
      check = InverseConsistentLandmarks::ReadIPLCerebellarPointTypes(lmrkfilename);
    }
    break;
    case UNKNOWN_LANDMARKS:
    {
      std::cout << "Unknown landmark file format!" << std::endl;
      return false;
    }
    default:
      itkGenericExceptionMacro(<< "Error. Invalid Landmark_File_Format type!");
  }
  return check;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadPointTypes(const std::string lmrkfilename,
                                                                           const int         XDim,
                                                                           const int         YDim,
                                                                           const int         ZDim,
                                                                           const int         TDim)
{
  Landmark_File_Format lmkff = QueryLandmarkFile(lmrkfilename);

  switch (lmkff)
  {
    case GEC_LANDMARKS:
    {
      InverseConsistentLandmarks::ReadGECPointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
    }
    break;
    case INTELLX_LANDMARKS:
    {
      InverseConsistentLandmarks::ReadIntellXPointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
    }
    break;
    case ANALYZE_LANDMARKS:
    {
      InverseConsistentLandmarks::ReadAnalyzePointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
    }
    break;
    case IPL_LANDMARKS:
    {
      InverseConsistentLandmarks::ReadIPLPointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
    }
    break;
    case IPL_TALAIRACH_BOUNDS_LANDMARKS:
    {
      InverseConsistentLandmarks::ReadIPLTalairachPointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
    }
    break;
    case IPL_CEREBELLAR_BOUNDS_LANDMARKS:
    {
      InverseConsistentLandmarks::ReadIPLCerebellarPointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
    }
    break;
    case UNKNOWN_LANDMARKS:
    {
      std::cout << "Unknown landmark file format!" << std::endl;
      return false;
    }
    default:
      itkGenericExceptionObject(<< "Error. Invalid Landmark_File_Format type!");
  }
  return true;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::WritePointTypes(const std::string          lmrkfilename,
                                                                            const Landmark_File_Format lmkff)
{
  bool writeresult;

  switch (lmkff)
  {
    case GEC_LANDMARKS:
    {
      writeresult = InverseConsistentLandmarks::WriteGECPointTypes(lmrkfilename);
    }
    break;
    case INTELLX_LANDMARKS:
    {
      writeresult = InverseConsistentLandmarks::WriteIntellXPointTypes(lmrkfilename);
    }
    break;
    case ANALYZE_LANDMARKS:
    {
      writeresult = InverseConsistentLandmarks::WriteAnalyzePointTypes(lmrkfilename);
    }
    break;
    case IPL_LANDMARKS:
    case IPL_TALAIRACH_BOUNDS_LANDMARKS:
    case IPL_CEREBELLAR_BOUNDS_LANDMARKS:
    {
      // writeresult = InverseConsistentLandmarks::WriteIPLPointTypes(
      // lmrkfilename );
      std::cout << "IPL_LANDMARKS NOT IMPLEMENTED FOR WRITING" << std::endl;
    }
    break;
    case UNKNOWN_LANDMARKS:
      std::cout << "Unknown output landmark file format!" << std::endl;
      writeresult = false;
  }
  return writeresult;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::PrintPointTypes(void) const
{
  itk::NumberToString<double> doubleConvert;

  printf("#IMAGEDIMS %6hu %6hu %6hu %6hu\n", ImageDims[0], ImageDims[1], ImageDims[2], ImageDims[3]);
  printf("#%11s %12s %12s %12s %12s %12s\n", "Name", "X", "Y", "Z", "T", "Weighting");
  for (typename InverseConsistentLandmarks::const_iterator mapiter = this->begin(); mapiter != this->end(); ++mapiter)
  {
    std::cout << mapiter->first;
    for (int space = 0; space < 12 - mapiter->first.size(); ++space)
    {
      std::cout << ' ';
    }
    std::cout << ' ';
    std::cout << doubleConvert(mapIter[0]) << ' ' << doubleConvert(mapIter[1]) << ' ' << doubleConvert(mapIter[2])
              << ' ' << mapiter->second.GetT() << ' ' << mapiter->second.GetWeighting() << std::endl;
    // printf( "%12s %12.6f %12.6f %12.6f %12.5f %12.6f\n", mapiter->first.c_str(),
    //         mapiter->second[0],
    //         mapiter->second[1],
    //         mapiter->second[2],
    //         mapiter->second.GetT(),
    //         mapiter->second.GetWeighting() );
  }
  if (this->size() > 0)
  {
    return true;
  }
  return false;
}

template <typename PointStorageType, typename PointSetType>
typename InverseConsistentLandmarks<PointStorageType, PointSetType>::Landmark_File_Format
InverseConsistentLandmarks<PointStorageType, PointSetType>::QueryLandmarkFile(const std::string lmrkfilename)
{
  /* These are the header formats for GEC, IntellX, and Analyze Landmark files
   * respectively
   * GEC Landmark File Format (typical extension .lmks)
   *
   *--------------------------------------------------------------------------------------
   * #GECLANDMARKS-v1.0
   * #IMAGEDIMS   XDIM   YDIM   ZDIM   TDIM
   * #Name            X            Y            Z            T
   *    Weighting
   * AC_Point 128 184 128 0 0
   * PC_Point 100 160 130 0 0
   * etc.
   *
   *--------------------------------------------------------------------------------------
   * IntellX Landmark File Format (typical extension .lmk)
   *
   *--------------------------------------------------------------------------------------
   * Landmarks-1.0
   * Number of points
   * "landmark identifier/description"
   * XPOINT YPOINT ZPOINT 1 1
   * "landmark identifier/description"
   * XPOINT YPOINT ZPOINT 1 1
   * etc.
   *
   *--------------------------------------------------------------------------------------
   * Analyze Landmark File Format (typical extension .log)
   *
   *--------------------------------------------------------------------------------------
   * #System Time
   * #Volume=name of volume
   * #
   * #
   * #Point
   * #X    Y    Z     Value       Label
   * #==== ==== ====  ==========  ==========
   *   129  277  128          14        Label
   *   150  100  164          32        Label
   * etc.
   *
   *--------------------------------------------------------------------------------------
   * For specializing which file format it is, the firstline #GECLANDMARKS-v1.0
   * will
   * specify a type GEC_LANDMARKS.  The tag Landmarks-1.0 will specify a
   * landmark
   * file of the type INTELLX_LANDMARKS.  The lines
   */
  // #Point
  // #X    Y    Z     Value       Label
  // will be used to specify a landmark file generated from an Analyze point log
  // session as the type ANALYZE_LANDMARKS.  Any other format which doesn't
  // match
  // will return UNKNOWN_LANDMARKS.

  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");

  if (tempfile == NULL)
  {
    std::cout << "ERROR: can not open " << lmrkfilename << std::endl;
    return UNKNOWN_LANDMARKS;
  }

  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return UNKNOWN_LANDMARKS;
  }

  if (strncmp(buffer, "#GECLANDMARKS-v1.0", 18) == 0)
  {
    return GEC_LANDMARKS;
  }

  if (strncmp(buffer, "Landmarks-1.0", 13) == 0)
  {
    return INTELLX_LANDMARKS;
  }

  if (strncmp(buffer, "IPL_HEADER_BEGIN", 16) == 0)
  {
    if (lmrkfilename.find("Talairach.bnd") != std::string::npos)
    {
      return IPL_TALAIRACH_BOUNDS_LANDMARKS;
    }
    else if (lmrkfilename.find("Cerebellum.bnd") != std::string::npos)
    {
      return IPL_CEREBELLAR_BOUNDS_LANDMARKS;
    }
    else
    {
      return IPL_LANDMARKS;
    }
  }
  // If it is not either of these two, try to read the next 4 lines
  // and check to see if it matches the Analyze point log file format
  for (int i = 0; i < 4; ++i)
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return UNKNOWN_LANDMARKS;
    }
  }

  if (strncmp(buffer, "#Point", 11) == 0)
  {
    // Try to match next line
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return UNKNOWN_LANDMARKS;
    }
    if (strncmp(buffer, "#X    Y    Z     Value       Label", 37) == 0)
    {
      return ANALYZE_LANDMARKS;
    }
  }
  return UNKNOWN_LANDMARKS;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadGECPointTypes(const std::string lmrkfilename)
{
  // std::cout << "Reading GEC Landmark File " << lmrkfilename << std::endl;

  this->erase(this->begin(), this->end());
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  constexpr unsigned short int LANDMARK_NAME_SIZE = 40;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char                         CurrentLandmarkName[LANDMARK_NAME_SIZE];
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  if (strncmp(buffer, "#GECLANDMARKS-v1.0", 18) != 0)
  {
    // Only files with a header of #GECLANDMARKS-v1.0 will be accepted
    return false;
  }
  // Read dimensions of image that landmarks are associated with
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  // Note CurrentLandmarkName is just a placeholder here!.
  // Note the hu is an indicator that this is a unsigned short int

  sscanf(buffer,
         "%s %hu %hu %hu %hu",
         CurrentLandmarkName,
         &(ImageDims[0]),
         &(ImageDims[1]),
         &(ImageDims[2]),
         &(ImageDims[3]));

  while (fgets(buffer, FILE_BUFFER_SIZE, tempfile) != NULL)
  {
    if (buffer[0] == '#') // This allows for comments to be added
                          // anywhere in the file
    {
      continue;
    }
    PointType TempPnt;

    double       val0;
    double       val1;
    double       val2;
    double       valT;
    double       valW;
    unsigned int numread =
      sscanf(buffer, "%s %lf %lf %lf %lf %lf", CurrentLandmarkName, &val0, &val1, &val2, &valT, &valW);
    if (numread < 6)
    {
      return false;
    }
    TempPnt[0] = static_cast<double>(val0);
    TempPnt[1] = static_cast<double>(val1);
    TempPnt[2] = static_cast<double>(val2);
    TempPnt.SetT(static_cast<double>(valT));
    TempPnt.SetWeighting(static_cast<double>(valW));

    if (ImageDims[0] <= TempPnt[0] || ImageDims[1] <= TempPnt[1] || ImageDims[2] <= TempPnt[2] ||
        ImageDims[3] <= TempPnt.GetT())
    {
      return false;
    }
    (*this)[CurrentLandmarkName] = TempPnt;
  }

  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadGECPointTypes(const std::string lmrkfilename,
                                                                              const int         XDim,
                                                                              const int         YDim,
                                                                              const int         ZDim,
                                                                              const int         TDim)
{
  const bool testread = this->InverseConsistentLandmarks::ReadGECPointTypes(lmrkfilename);

  this->InverseConsistentLandmarks::rescale(XDim, YDim, ZDim, TDim);
  return testread;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIntellXPointTypes(const std::string lmrkfilename)
{
  //  std::cout << "Reading IntellX Landmark File " << lmrkfilename <<
  // std::endl;
  this->erase(this->begin(), this->end());
  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 128;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  std::string                  CurrentLandmarkName;
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  if (strncmp(buffer, "Landmarks-1.0", 13) != 0)
  {
    // Only files with a header of Landmarks-1.0 will be accepted
    return false;
  }
  // Read Number of Landmark points
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  const int NumberOfPoints = std::stoi(buffer);
  for (int p = 0; p < NumberOfPoints; ++p)
  {
    // Extract Landmark Name
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
    CurrentLandmarkName = std::string(buffer);
    // CurrentLandmarkName.replace(CurrentLandmarkName.find("\""),1,"");
    // CurrentLandmarkName.replace(CurrentLandmarkName.find("\""),1,"");
    std::string NewLandmarkName;
    for (unsigned int iter = 0; iter < CurrentLandmarkName.length(); ++iter)
    {
      // std::cout << "CurrentLandmarkName[" << iter << "] = " <<
      // CurrentLandmarkName[iter] << std::endl;
      if (CurrentLandmarkName[iter] == '\n' || CurrentLandmarkName[iter] == '\r' || CurrentLandmarkName[iter] == '\t' ||
          CurrentLandmarkName[iter] == '\"')
      {
        continue;
      }
      else if (CurrentLandmarkName[iter] == ' ')
      {
        NewLandmarkName += "_";
      }
      else
      {
        NewLandmarkName += CurrentLandmarkName[iter];
      }
      // std::cout << "NewLandmarkName = " << NewLandmarkName << std::endl;
    }

    // Extract Coordinates
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
    PointType TempPnt;
    {
      unsigned int numread = sscanf(buffer, "%lf %lf %lf", &TempPnt[0], &TempPnt[1], &TempPnt[2]);
      if (numread < 3)
      {
        return false;
      }
    }
    // Convert to Analyze coordinate system...Note the -1 is needed for correct
    // coord
    // Ex. 0,0,0 in Warp = 127,127,127 in a 128x128x128 analyze image.
    // Assumes the landmark dimensions have been set before used in this
    // function
    TempPnt[0] = static_cast<PointStorageType>(ImageDims[0]) - TempPnt[0];
    TempPnt[1] = static_cast<PointStorageType>(ImageDims[1]) - TempPnt[1];
    TempPnt[2] = static_cast<PointStorageType>(ImageDims[2]) - TempPnt[2];
    // DDEBUG, this is set to one when dimension is one, which goofs when you
    // write out
    // a lmk file as a ilmks and try to read the ilmks back in.  Since no t
    // coordinates
    // will be used in the short run, this will be implemented
    // TempPnt.setT() = static_cast<PointStorageType>( ImageDims[3] ) -
    // TempPnt[3];
    TempPnt.SetT(0.0F);

    // Set Weighting to default value of 1.0F;
    TempPnt.SetWeighting(1.0F);
    (*this)[NewLandmarkName] = TempPnt;
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIntellXPointTypes(const std::string lmrkfilename,
                                                                                  const int         XDim,
                                                                                  const int         YDim,
                                                                                  const int         ZDim,
                                                                                  const int         TDim)
{
  // std::cout << "Reading IntellX Landmark File " << lmrkfilename << std::endl;
  // std::cout << "Specifying dimensions " << XDim << ", " << YDim << ", " <<
  // ZDim << ", " << TDim << "." << std::endl;

  if (XDim < 0)
  {
    std::cout << "Error. X dimension " << XDim << " given is less than 0" << std::endl;
    return false;
  }
  if (YDim < 0)
  {
    std::cout << "Error. Y dimension " << YDim << " given is less than 0" << std::endl;
    return false;
  }
  if (ZDim < 0)
  {
    std::cout << "Error. Z dimension " << ZDim << " given is less than 0" << std::endl;
    return false;
  }
  if (TDim < 0)
  {
    std::cout << "Error. T dimension " << TDim << " given is less than 0" << std::endl;
    return false;
  }

  // Delete all old landmarks
  this->erase(this->begin(), this->end());
  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 128;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  std::string                  CurrentLandmarkName;
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  if (strncmp(buffer, "Landmarks-1.0", 13) != 0)
  {
    // Only files with a header of Landmarks-1.0 will be accepted
    return false;
  }
  // Read Number of Landmark points
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  const int NumberOfPoints = std::stoi(buffer);

  // Assigning the dimensions from the parameters
  ImageDims[0] = XDim;
  ImageDims[1] = YDim;
  ImageDims[2] = ZDim;
  ImageDims[3] = TDim;
  for (int p = 0; p < NumberOfPoints; ++p)
  {
    // Extract Landmark Name
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }

    CurrentLandmarkName = std::string(buffer);

    std::string NewLandmarkName;
    for (unsigned int iter = 0; iter < CurrentLandmarkName.length(); ++iter)
    {
      // std::cout << "CurrentLandmarkName[" << iter << "] = " <<
      // CurrentLandmarkName[iter] << std::endl;
      if (CurrentLandmarkName[iter] == '\n' || CurrentLandmarkName[iter] == '\r' || CurrentLandmarkName[iter] == '\t' ||
          CurrentLandmarkName[iter] == '\"')
      {
        continue;
      }
      else if (CurrentLandmarkName[iter] == ' ')
      {
        NewLandmarkName += "_";
      }
      else
      {
        NewLandmarkName += CurrentLandmarkName[iter];
      }
    }

    // Extract Coordinates
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
    PointType TempPnt;
    {
      double       val0;
      double       val1;
      double       val2;
      unsigned int numread = sscanf(buffer, "%lf %lf %lf", &val0, &val1, &val2);
      if (numread < 3)
      {
        return false;
      }
      TempPnt[0] = static_cast<double>(val0);
      TempPnt[1] = static_cast<double>(val1);
      TempPnt[2] = static_cast<double>(val2);
    }
    // Convert to Analyze coordinate system...Note the -1 is needed for correct
    // coord
    // Ex. 0,0,0 in Warp = 127,127,127 in a 128x128x128 analyze image.
    TempPnt[0] = static_cast<PointStorageType>(ImageDims[0]) - TempPnt[0];
    TempPnt[1] = static_cast<PointStorageType>(ImageDims[1]) - TempPnt[1];
    TempPnt[2] = static_cast<PointStorageType>(ImageDims[2]) - TempPnt[2];

    TempPnt.SetT(0.0F);

    // Set Weighting to default value of 1.0F;
    TempPnt.SetWeighting(1.0F);
    (*this)[NewLandmarkName] = TempPnt;
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadAnalyzePointTypes(const std::string lmrkfilename)
{
  std::cout << "Reading Analyze Landmark File " << lmrkfilename << std::endl;

  // Delete all old landmarks
  this->erase(this->begin(), this->end());

  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  constexpr unsigned short int LANDMARK_NAME_SIZE = 40;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char                         CurrentLandmarkName[LANDMARK_NAME_SIZE];
  char *                       status = 0;
  // Read Header Information
  for (int i = 0; i < 5; ++i)
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
  }

  if (strncmp(buffer, "#Point", 11) == 0)
  {
    // Try to match next line
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
    if (strncmp(buffer, "#X    Y    Z     Value       Label", 37) != 0)
    {
      return false;
    }
  }
  else
  {
    return false;
  }

  // Reading the line with " ====  ====  ==== ..etc"
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }

  int mycount = 1;
  while (fgets(buffer, FILE_BUFFER_SIZE, tempfile) != NULL)
  {
    int       grayvalue;
    PointType TempPnt;
    {
      double val0;
      double val1;
      double val2;
      int    numread = sscanf(buffer, "%lf %lf %lf %d %s", &val0, &val1, &val2, &grayvalue, CurrentLandmarkName);
      TempPnt[0] = static_cast<double>(val0);
      TempPnt[1] = static_cast<double>(val1);
      TempPnt[2] = static_cast<double>(val2);
      if (numread == -1)
      {
        continue;
      }

      if (numread < 4) // Need at least an X,Y,Z coordinate
      {
        return false;
      }
      if (numread == 4) // No label was found
      {
        memset(CurrentLandmarkName, 0, LANDMARK_NAME_SIZE);
      }
    }

    if (ImageDims[0] <= TempPnt[0])
    {
      return false;
    }
    if (ImageDims[1] <= TempPnt[1])
    {
      return false;
    }
    if (ImageDims[2] <= TempPnt[2])
    {
      return false;
    }

    // Here is how the naming definition is defined, if the label is not found,
    // the
    // landmark identifier is Landmark_N where N is the current number of
    // landmarks
    // read in.  If the landmark label is found in the current map, then the
    // label is renamed with the old label plus a "_1".  If that tag exists,
    // then
    // the next label is the old label with "_2".  This continues until a
    // unique name is found for the landmark.

    std::string CurrentName(CurrentLandmarkName);

    if (CurrentName.size() == 0)
    {
      char numberbuffer[16];
      sprintf(numberbuffer, "%d", mycount);
      CurrentName = "Landmark_" + std::string(numberbuffer);
    }

    // Checking to make sure no elements in the map of landmark points are
    // duplicated.
    {
      int         label = 1;
      std::string OriginalLabel(CurrentName);

      char labelbuffer[16];
      while (this->find(CurrentName) != this->end())
      {
        sprintf(labelbuffer, "%d", label);
        CurrentName = OriginalLabel + "_" + std::string(labelbuffer);
        ++label;
      }
    }
    TempPnt.SetWeighting(1.0F);
    (*this)[CurrentName] = TempPnt;
    ++mycount;
  }

  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadAnalyzePointTypes(const std::string lmrkfilename,
                                                                                  const int         XDim,
                                                                                  const int         YDim,
                                                                                  const int         ZDim,
                                                                                  const int         TDim)
{
  std::cout << "Reading Analyze Landmark File " << lmrkfilename << std::endl;
  std::cout << "Specifying dimensions " << XDim << ", " << YDim << ", " << ZDim << ", " << TDim << "." << std::endl;

  if (XDim < 0)
  {
    std::cout << "Error. X dimension " << XDim << " given is less than 0" << std::endl;
    return false;
  }
  if (YDim < 0)
  {
    std::cout << "Error. Y dimension " << YDim << " given is less than 0" << std::endl;
    return false;
  }
  if (ZDim < 0)
  {
    std::cout << "Error. Z dimension " << ZDim << " given is less than 0" << std::endl;
    return false;
  }
  if (TDim < 0)
  {
    std::cout << "Error. T dimension " << TDim << " given is less than 0" << std::endl;
    return false;
  }

  // Delete all old landmarks
  this->erase(this->begin(), this->end());

  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  constexpr unsigned short int LANDMARK_NAME_SIZE = 40;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char                         CurrentLandmarkName[LANDMARK_NAME_SIZE];
  char *                       status = 0;
  // Read Header Information
  for (int i = 0; i < 5; ++i)
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
  }

  if (strncmp(buffer, "#Point", 11) == 0)
  {
    // Try to match next line
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      return false;
    }
    if (strncmp(buffer, "#X    Y    Z     Value       Label", 37) != 0)
    {
      return false;
    }
  }
  else
  {
    return false;
  }

  // Reading the line with " ====  ====  ==== ..etc"
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }

  // Using the function to assign the landmark image dimensions
  ImageDims[0] = XDim;
  ImageDims[1] = YDim;
  ImageDims[2] = ZDim;
  ImageDims[3] = TDim;

  int count = 1;
  while (fgets(buffer, FILE_BUFFER_SIZE, tempfile) != NULL)
  {
    int       grayvalue;
    PointType TempPnt;
    {
      double val0;
      double val1;
      double val2;
      int    numread = sscanf(buffer, "%lf %lf %lf %d %s", &val0, &val1, &val2, &grayvalue, CurrentLandmarkName);
      TempPnt[0] = static_cast<double>(val0);
      TempPnt[1] = static_cast<double>(val1);
      TempPnt[2] = static_cast<double>(val2);

      if (numread == -1)
      {
        // No values read in from sscanf so continue
        continue;
      }

      if (numread < 4) // Need at least an X,Y,Z coordinate
      {
        return false;
      }

      if (numread == 4) // No label was found
      {
        memset(CurrentLandmarkName, 0, LANDMARK_NAME_SIZE);
      }

      // std::cout << "Number read = " << numread << std::endl;
      // std::cout << "CurrentLandmarkName = " << CurrentLandmarkName <<
      // std::endl;
    }

    if (ImageDims[0] <= TempPnt[0])
    {
      return false;
    }
    if (ImageDims[1] <= TempPnt[1])
    {
      return false;
    }
    if (ImageDims[2] <= TempPnt[2])
    {
      return false;
    }

    // Here is how the naming definition is defined, if the label is not found,
    // the
    // landmark identifier is Landmark_N where N is the current number of
    // landmarks
    // read in.  If the landmark label is found in the current map, then the
    // label is renamed with the old label plus a "_1".  If that tag exists,
    // then
    // the next label is the old label with "_2".  This continues until a
    // unique name is found for the landmark.

    std::string CurrentName(CurrentLandmarkName);

    if (CurrentName.size() == 0)
    {
      char numberbuffer[16];
      sprintf(numberbuffer, "%d", count);
      CurrentName = "Landmark_" + std::string(numberbuffer);
      // std::cout << "Landmark label = " << CurrentName << std::endl;
    }

    // Checking to make sure no elements in the map of landmark points are
    // duplicated.
    {
      int         label = 1;
      std::string OriginalLabel(CurrentName);

      char labelbuffer[16];
      while (this->find(CurrentName) != this->end())
      {
        sprintf(labelbuffer, "%d", label);
        CurrentName = OriginalLabel + "_" + std::string(labelbuffer);
        ++label;
      }

      // std::cout << "CurrentName = " << CurrentName << std::endl;
    }
    TempPnt.SetWeighting() = 1.0F;
    (*this)[CurrentName] = TempPnt;
    ++count;
  }

  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::WriteGECPointTypes(const std::string lmrkfilename) const
{
  // Open the files
  FILE *                      tempfile = fopen(lmrkfilename.c_str(), "w");
  itk::NumberToString<double> doubleConvert;

  if (tempfile == NULL)
  {
    return false;
  }
  fprintf(tempfile, "#GECLANDMARKS-v1.0\n");
  fprintf(tempfile, "#IMAGEDIMS %6hu %6hu %6hu %6hu\n", ImageDims[0], ImageDims[1], ImageDims[2], ImageDims[3]);
  // Note: This is just a comment
  fprintf(tempfile, "#%11s %12s %12s %12s %12s %12s\n", "Name", "X", "Y", "Z", "T", "Weighting");
  for (typename InverseConsistentLandmarks::const_iterator mapiter = this->begin(); mapiter != this->end(); ++mapiter)
  {
    //    fprintf( tempfile, "%12s %12.6f %12.6f %12.6f %12.5f %12.6f\n",
    // mapiter->first.c_str(),
    // mapiter->second[0],
    // mapiter->second[1],
    // mapiter->second[2],
    // mapiter->second.GetT(),
    // mapiter->second.GetWeighting() );
    std::string s = mapiter->first;
    s += ' ';
    for (unsigned dir = 0; dir < 3; ++dir)
    {
      s += doubleConvert(mapiter->second[dir]);
      s += ' ';
    }
    s += doubleConvert(mapiter->second.GetT());
    s += ' ';
    s += doubleConvert(mapIter->GetWeighting);
    s += '\n';
    fputs(s.c_str(), tempFile);
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::WriteIntellXPointTypes(const std::string lmrkfilename) const
{
  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "w");

  if (tempfile == NULL)
  {
    return false;
  }
  fprintf(tempfile, "Landmarks-1.0\n");
  fprintf(tempfile, "%d\n", static_cast<int>(this->size()));
  for (typename InverseConsistentLandmarks::const_iterator mapiter = this->begin(); mapiter != this->end(); ++mapiter)
  {
    std::string temp("\"");

    temp += mapiter->first;
    temp += "\"";
    // Notice the conversion back to IntellX coordinate system.
    // Note: -1 is needed here to get correct coordinate.
    fprintf(tempfile,
            "%s\n%f\t%f\t%f\t1\t1\n",
            temp.c_str(),
            static_cast<PointStorageType>(ImageDims[0]) - mapiter->second[0],
            static_cast<PointStorageType>(ImageDims[1]) - mapiter->second[1],
            static_cast<PointStorageType>(ImageDims[2]) - mapiter->second[2]);
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::WriteAnalyzePointTypes(const std::string lmrkfilename) const
{
  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "w");

  if (tempfile == NULL)
  {
    return false;
  }

  // Getting the current time information
  time_t      curtime;
  struct tm * loctime;
  curtime = time(NULL);
  loctime = localtime(&curtime);
  std::string date(asctime(loctime));

  date = std::string("#") + date;
  fprintf(tempfile, date.c_str());
  fprintf(tempfile, "#Volume=\n");
  fprintf(tempfile, "#\n");
  fprintf(tempfile, "#\n");
  fprintf(tempfile, "#Point\n");
  fprintf(tempfile, "#X    Y    Z     Value       Label\n");
  fprintf(tempfile, "#==== ==== ====  ==========  ==========\n");
  for (typename InverseConsistentLandmarks::const_iterator mapiter = this->begin(); mapiter != this->end(); ++mapiter)
  {
    // Graylevel value is assumed to not be necessary so it's set to 0
    fprintf(tempfile,
            "%5d%5d%5d%12d    %s\n",
            static_cast<int>(mapiter->second[0]),
            static_cast<int>(mapiter->second[1]),
            static_cast<int>(mapiter->second[2]),
            0,
            mapiter->first.c_str());
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ConcatLandmarks(InverseConsistentLandmarks & newlmks)
{
  // Copy landmarks from newlmks to the current landmark file, but
  // ensuring that the keys are not identical

  if (this == &newlmks)
  {
    std::cout << "Error. Can not concatenate a landmark object with itself!" << std::endl;
    return false;
  }

  // Check to make sure the dimensions match
  bool sizematch = true;
  if (this->getXDim() != newlmks.getXDim())
  {
    std::cout << "Error. X dimensions of landmark objects do not match in ConcatLandmarks!" << std::endl;
    sizematch = false;
  }
  if (this->getYDim() != newlmks.getYDim())
  {
    std::cout << "Error. Y dimensions of landmark objects do not match in ConcatLandmarks!" << std::endl;
    sizematch = false;
  }
  if (this->getZDim() != newlmks.getZDim())
  {
    std::cout << "Error. Z dimensions of landmark objects do not match in ConcatLandmarks!" << std::endl;
    sizematch = false;
  }
  if (this->getTDim() != newlmks.getTDim())
  {
    std::cout << "Error. T dimensions of landmark objects do not match in ConcatLandmarks!" << std::endl;
    sizematch = false;
  }
  if (sizematch == false)
  {
    return sizematch;
  }

  // This will store all the labels
  std::map<std::string, int> labelcount;
  std::string                CurrentName;
  char                       labelbuffer[16];

  if (newlmks.size() == 0)
  {
    return false;
  }
  // Initialize the map with the labels found in this landmark object
  for (typename InverseConsistentLandmarks::iterator iter = this->begin(); iter != this->end(); ++iter)
  {
    labelcount[iter->first] = 0;
  }
  for (typename InverseConsistentLandmarks::iterator iter = newlmks.begin(); iter != newlmks.end(); ++iter)
  {
    if (labelcount.find(iter->first) != labelcount.end())
    {
      // Increment the counter and apply that name with the new tag
      labelcount[iter->first] = labelcount[iter->first] + 1;
      sprintf(labelbuffer, "%d", labelcount.find(iter->first)->second);
      CurrentName = iter->first + "_" + std::string(labelbuffer);
    }
    else
    {
      labelcount[iter->first] = 0;
      CurrentName = iter->first;
    }
    // std::cout << "CurrentName = " << CurrentName << std::endl;
    (*this)[CurrentName] = iter->second;
  }

  return true;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ConcatLandmarks(const std::string lmrkfilename)
{
  // Copy landmarks from newlmks to the current landmark file
  InverseConsistentLandmarks temp;

  temp.ReadPointTypes(lmrkfilename);
  return this->ConcatLandmarks(temp);
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ConcatLandmarks(const std::string lmrkfilename,
                                                                            const int         XDim,
                                                                            const int         YDim,
                                                                            const int         ZDim,
                                                                            const int         TDim)
{
  // Copy landmarks from newlmks to the current landmark file
  InverseConsistentLandmarks temp;

  temp.ReadPointTypes(lmrkfilename, XDim, YDim, ZDim, TDim);
  return this->ConcatLandmarks(temp);
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::RemoveClosePoints(const PointStorageType distance)
{
  bool removedpoint = false;

  typename InverseConsistentLandmarks::iterator iter;
  typename InverseConsistentLandmarks::iterator jter;
  for (iter = this->begin(); iter != this->end(); ++iter)
  {
    // std::cout << "iter->first = " << iter->first << std::endl;
    typename InverseConsistentLandmarks::iterator ipp = iter;
    ++ipp;
    jter = ipp;

    while (jter != this->end())
    {
      // std::cout << "jter->first = " << jter->first << std::endl;
      PointStorageType ix = iter->second[0];
      PointStorageType iy = iter->second[1];
      PointStorageType iz = iter->second[2];
      PointStorageType jx = jter->second[0];
      PointStorageType jy = jter->second[1];
      PointStorageType jz = jter->second[2];

      PointStorageType value = sqrt((jx - ix) * (jx - ix) + (jy - iy) * (jy - iy) + (jz - iz) * (jz - iz));
      if (value < distance)
      {
        // You lose the iterator when you issue an erase, so saving jter
        typename InverseConsistentLandmarks::iterator dter = jter;
        removedpoint = true;
        this->erase(dter->first);

        // The iterator has stayed the same, but the list has moved back one, so
        // an
        // extra increment is needed.
        ++jter;
      }
      ++jter;
    }
  }

  return removedpoint;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::RemoveUnmatchedPoints(InverseConsistentLandmarks & tempmap1)
{
  for (typename InverseConsistentLandmarks::iterator mapiter = this->begin(); mapiter != this->end(); ++mapiter)
  {
    typename InverseConsistentLandmarks::iterator mapiter1 = tempmap1.find(mapiter->first);
    // Item NOT Found
    if (mapiter1 == tempmap1.end())
    {
      this->erase(mapiter->first);
    }
  }

  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::rescale(const int newx,
                                                                    const int newy,
                                                                    const int newz,
                                                                    const int newt)
{
  assert(newx > 0);
  assert(newy > 0);
  assert(newz > 0);
  assert(newt > 0);

  const PointStorageType scalex = static_cast<PointStorageType>(newx) / static_cast<PointStorageType>(ImageDims[0]);
  const PointStorageType scaley = static_cast<PointStorageType>(newy) / static_cast<PointStorageType>(ImageDims[1]);
  const PointStorageType scalez = static_cast<PointStorageType>(newz) / static_cast<PointStorageType>(ImageDims[2]);
  const PointStorageType scalet = static_cast<PointStorageType>(newt) / static_cast<PointStorageType>(ImageDims[3]);

  ImageDims[0] = static_cast<unsigned short>(newx);
  ImageDims[1] = static_cast<unsigned short>(newy);
  ImageDims[2] = static_cast<unsigned short>(newz);
  ImageDims[3] = static_cast<unsigned short>(newt);
  std::cout << "Dims: " << ImageDims[0] << " " << ImageDims[1] << " " << ImageDims[2] << " " << ImageDims[3] << "\n";
  std::cout << "scale: " << scalex << " " << scaley << " " << scalez << " " << scalet << "\n";
  for (typename InverseConsistentLandmarks::iterator mapiter = this->begin(); mapiter != this->end(); ++mapiter)
  {
    mapiter->second[0] = mapiter->second[0] * scalex;
    mapiter->second[1] = mapiter->second[1] * scaley;
    mapiter->second[2] = mapiter->second[2] * scalez;
    mapiter->second.SetT(mapiter->second.GetT() * scalet);
  }

  return true;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIPLPointTypes(const std::string lmrkfilename)
{
  std::cout << "Reading IPL Landmark File " << lmrkfilename << std::endl;

  // Delete all old landmarks
  this->erase(this->begin(), this->end());
  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  constexpr unsigned short int LANDMARK_NAME_SIZE = 40;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char                         CurrentLandmarkName[LANDMARK_NAME_SIZE];
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  if (strncmp(buffer, "IPL_HEADER_BEGIN", 16) != 0)
  {
    // Only files with a header of IPL_HEADER_BEGIN will be accepted
    return false;
  }
  // Read dimensions of image that landmarks are associated with
  // skip a head to landmark definitions
  status = buffer;
  while ((strncmp(buffer, "LANDMARK_HEADER_BEGIN", 21) != 0) && (status != NULL))
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before LANDMARK_HEADER_BEGIN found" << std::endl;
    }
  }

  bool xisset = false;
  bool yisset = false;
  bool zisset = false;
  ImageRes[0] = 1.0;
  ImageRes[1] = 1.0;
  ImageRes[2] = 1.0;
  ImageRes[3] = 1.0;

  // process looking for dims
  status = buffer;
  while ((strncmp(buffer, "LANDMARK_HEADER_END", 19) != 0) && (status != NULL))
  {
    // Look for x
    if (strncmp(buffer, "LANDMARK_X_SIZE", 15) == 0)
    {
      sscanf(buffer, "%s %hu", CurrentLandmarkName, &(ImageDims[0]));
      xisset = true;
    }
    // Look for y
    if (strncmp(buffer, "LANDMARK_Y_SIZE", 15) == 0)
    {
      sscanf(buffer, "%s %hu", CurrentLandmarkName, &(ImageDims[1]));
      yisset = true;
    }
    // Look for z
    if (strncmp(buffer, "LANDMARK_Z_SIZE", 15) == 0)
    {
      sscanf(buffer, "%s %hu", CurrentLandmarkName, &(ImageDims[2]));
      zisset = true;
    }
    if (strncmp(buffer, "LANDMARK_X_RESOLUTION", 20) == 0)
    {
      sscanf(buffer, "%s %f", CurrentLandmarkName, &(ImageRes[0]));
    }
    if (strncmp(buffer, "LANDMARK_Y_RESOLUTION", 20) == 0)
    {
      sscanf(buffer, "%s %f", CurrentLandmarkName, &(ImageRes[1]));
    }
    if (strncmp(buffer, "LANDMARK_Z_RESOLUTION", 20) == 0)
    {
      sscanf(buffer, "%s %f", CurrentLandmarkName, &(ImageRes[2]));
    }
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before LANDMARK_HEADER_END found" << std::endl;
    }
  }

  ImageDims[3] = 1;
  if (!(xisset && yisset && zisset))
  {
    std::cout << "ERROR must set all LANDMARK_[XYZ]_SIZE values" << std::endl;
  }

  // Skip ahead to data
  status = buffer;
  while ((strncmp(buffer, "IPL_HEADER_END", 14) != 0) && (status != NULL))
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before IPL_HEADER_END found" << std::endl;
    }
  }

  // Actually read the data
  while (fgets(buffer, FILE_BUFFER_SIZE, tempfile) != NULL)
  {
    if (buffer[0] == '#') // This allows for comments to be added
                          // anywhere in the file
    {
      continue;
    }
    PointType TempPnt;
    {
      // DO I NEED TO INVERT Y HERE? SMP
      // const double invertY = static_cast<double> (ImageDims[1]-1);
      unsigned int numread =
        sscanf(buffer, "%s %lf %lf %lf ", CurrentLandmarkName, &TempPnt[0], &TempPnt[1], &TempPnt[2]);
      if (numread != 4)
      {
        std::cout << "ERROR invalid number of args read at landmark " << CurrentLandmarkName << std::endl;
        return false;
      }
    }
    TempPnt.SetT(0.0F);
    TempPnt.SetWeighting(1.0F);
    if ((ImageDims[0] <= TempPnt[0] * ImageRes[0]) || (ImageDims[1] <= TempPnt[1] * ImageRes[1]) ||
        (ImageDims[2] <= TempPnt[2] * ImageRes[2]) || (ImageDims[3] <= TempPnt.GetT() * ImageRes[3]))
    {
      std::cout << "Point Set Outside of Image Dimensions." << std::endl;
      return false;
    }
    (*this)[CurrentLandmarkName] = TempPnt;
  }

  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIPLPointTypes(const std::string lmrkfilename,
                                                                              const int         XDim,
                                                                              const int         YDim,
                                                                              const int         ZDim,
                                                                              const int         TDim)
{
  const bool testread = this->InverseConsistentLandmarks::ReadIPLPointTypes(lmrkfilename);

  this->InverseConsistentLandmarks::rescale(XDim, YDim, ZDim, TDim);
  return testread;
}

// ==================Talairach Bounds to Landmark Points===================
template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::process_bnd_point(const std::string &  CurrentLandmarkName,
                                                                              const char *         buffer,
                                                                              const unsigned short local_ImageDims[4],
                                                                              const float /*local_ImageRes */[4],
                                                                              PointType & ModifiedPoint)
{
  PointStorageType tmp[3];

  {
    double       val0;
    double       val1;
    double       val2;
    unsigned int numread = sscanf(buffer, "%lf %lf %lf ", &val0, &val1, &val2);
    if (numread != 3)
    {
      std::cout << "ERROR invalid number of args read at landmark " << CurrentLandmarkName << std::endl;
      return false;
    }
    tmp[0] = val0;
    tmp[1] = val1;
    tmp[2] = val2;
  }

  for (int i = 0; i < 3; ++i)
  {
    ModifiedPoint[i] = tmp[i];
  }

  // Invert Y Readings for our coordinate system.
  ModifiedPoint[1] = local_ImageDims[1] - 1 - ModifiedPoint[1];
  ModifiedPoint.SetT(0.0F);
  ModifiedPoint.SetWeighting(1.0F);
  if (local_ImageDims[0] <= ModifiedPoint[0])
  {
    return false;
  }
  if (local_ImageDims[1] <= ModifiedPoint[1])
  {
    return false;
  }
  if (local_ImageDims[2] <= ModifiedPoint[2])
  {
    return false;
  }
  if (local_ImageDims[3] <= ModifiedPoint.GetT())
  {
    return false;
  }
  for (int i = 0; i < 3; ++i)
  {
    // ModifiedPoint[i]*=local_ImageRes[i];
  }
  return true;
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIPLTalairachPointTypes(const std::string lmrkfilename)
{
  // std::cout << "Reading IPL Talairach Bounds as  Landmark File " <<
  // lmrkfilename << std::endl;

  // Delete all old landmarks
  this->erase(this->begin(), this->end());
  // Open the files
  std::cout << "  File: " << lmrkfilename.c_str() << std::endl;
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");
  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  if (strncmp(buffer, "IPL_HEADER_BEGIN", 16) != 0)
  {
    // Only files with a header of IPL_HEADER_BEGIN will be accepted
    return false;
  }
  // Read dimensions of image that landmarks are associated with
  // skip a head to landmark definitions
  status = buffer;
  while ((strncmp(buffer, "TALAIRACH_PARAMETER_HEADER_BEGIN", 32) != 0) && (status != NULL))
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    // std::cout << "LINE: " << buffer << std::endl;
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_BEGIN found" << std::endl;
    }
  }

  bool xisset = false;
  bool yisset = false;
  bool zisset = false;
  // process looking for dims
  status = buffer;
  while ((strncmp(buffer, "TALAIRACH_PARAMETER_HEADER_END", 30) != 0) && (status != NULL))
  {
    constexpr unsigned short int TALAIRACH_PARAMETER_NAME_SIZE = 40;
    char                         Dummy[TALAIRACH_PARAMETER_NAME_SIZE];
    // Look for x
    if (strncmp(buffer, "TALAIRACH_PARAMETER_X_SIZE", 26) == 0)
    {
      sscanf(buffer, "%s %hu", Dummy, &(ImageDims[0]));
      xisset = true;
    }
    // Look for y
    if (strncmp(buffer, "TALAIRACH_PARAMETER_Y_SIZE", 26) == 0)
    {
      sscanf(buffer, "%s %hu", Dummy, &(ImageDims[1]));
      yisset = true;
    }
    // Look for z
    if (strncmp(buffer, "TALAIRACH_PARAMETER_Z_SIZE", 26) == 0)
    {
      sscanf(buffer, "%s %hu", Dummy, &(ImageDims[2]));
      zisset = true;
    }
    if (strncmp(buffer, "TALAIRACH_PARAMETER_X_RESOLUTION", 26) == 0)
    {
      sscanf(buffer, "%s %f", Dummy, &(ImageRes[0]));
    }
    // Look for y
    if (strncmp(buffer, "TALAIRACH_PARAMETER_Y_RESOLUTION", 26) == 0)
    {
      sscanf(buffer, "%s %f", Dummy, &(ImageRes[1]));
    }
    // Look for z
    if (strncmp(buffer, "TALAIRACH_PARAMETER_Z_RESOLUTION", 26) == 0)
    {
      sscanf(buffer, "%s %f", Dummy, &(ImageRes[2]));
    }
    // if(!(xisset && yisset && zisset))
    // {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    // }
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
  }

  ImageDims[3] = 1;
  ImageRes[3] = 1.0;
  if (!(xisset && yisset && zisset))
  {
    std::cout << "ERROR must set all TALAIRACH_PARAMETER_[XYZ]_SIZE values" << std::endl;
  }

  // Skip ahead to data
  status = buffer;
  while ((strncmp(buffer, "IPL_HEADER_END", 14) != 0) && (status != NULL))
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before IPL_HEADER_END found" << std::endl;
    }
  }

  // Actually read the data
  // First read the SLA point
  {
    PointType AC_Point; // Antieror Commisure
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("AC_Point", buffer, ImageDims, ImageRes, AC_Point) == false)
    {
      return false;
    }

    PointType PC_Point; // Posterior Commisure
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("PC_Point", buffer, ImageDims, ImageRes, PC_Point) == false)
    {
      return false;
    }

    PointType SLAPoint; // Superior Left Anterior
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("SLAPoint", buffer, ImageDims, ImageRes, SLAPoint) == false)
    {
      return false;
    }

    PointType IRPPoint; // Inferior Right Posterior
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("IRPPoint", buffer, ImageDims, ImageRes, IRPPoint) == false)
    {
      return false;
    }

    // Process to extend 4 points to all possible combinations
    PointType TempPnt = SLAPoint;
    {
      std::cout << "  Coordinates:: SLA " << SLAPoint << "   IRP " << IRPPoint << "   PC " << PC_Point << "   AC "
                << AC_Point << std::endl;
      std::string CurrentLandmarkName = "AC-Point";
      TempPnt[0] = AC_Point[0];
      TempPnt[1] = AC_Point[1];
      TempPnt[2] = AC_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "PC-Point";
      TempPnt[0] = PC_Point[0];
      TempPnt[1] = PC_Point[1];
      TempPnt[2] = PC_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;

      /* Inferior Points */
      CurrentLandmarkName = "IRP";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "IRA";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "ILP";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "ILA";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;

      /* Superior Points */
      CurrentLandmarkName = "SLA";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "SLP";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "SRA";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "SRP";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;

      /*
       *  CurrentLandmarkName="Tal_SLAPoint";    TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SLPPoint";    TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SRAPoint";    TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SRPPoint";    TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_S_PC_APoint"; TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_S_PC_PPoint"; TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SL_AC_Point"; TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SR_AC_Point"; TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SL_PC_Point"; TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_SR_PC_Point"; TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_S_AC_Point";  TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_S_PC_Point";  TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=SLAPoint[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *
       *  CurrentLandmarkName="Tal_PC_LAPoint";    TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_LPPoint";    TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_RAPoint";    TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_RPPoint";    TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_PC_APoint"; TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_PC_PPoint"; TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_L_AC_Point"; TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_R_AC_Point"; TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_L_PC_Point"; TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_R_PC_Point"; TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_AC_Point";      TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=AC_Point[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_PC_Point";      TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=PC_Point[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *
       *  CurrentLandmarkName="Tal_ILAPoint";    TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_ILPPoint";    TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_IRAPoint";    TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_IRPPoint";    TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_I_PC_APoint"; TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=SLAPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_I_PC_PPoint"; TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=IRPPoint[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_IL_AC_Point"; TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_IR_AC_Point"; TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_IL_PC_Point"; TempPnt[0]=SLAPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_IR_PC_Point"; TempPnt[0]=IRPPoint[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_I_AC_Point";  TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=AC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       *  CurrentLandmarkName="Tal_I_PC_Point";  TempPnt[0]=PC_Point[0];
       *  TempPnt[1]=IRPPoint[1]; TempPnt[2]=PC_Point[2];
       *     (*this)[CurrentLandmarkName]=TempPnt;
       */
    }
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIPLTalairachPointTypes(const std::string lmrkfilename,
                                                                                       const int         XDim,
                                                                                       const int         YDim,
                                                                                       const int         ZDim,
                                                                                       const int         TDim)
{
  const bool testread = this->InverseConsistentLandmarks::ReadIPLPointTypes(lmrkfilename);

  this->InverseConsistentLandmarks::rescale(XDim, YDim, ZDim, TDim);
  return testread;
}

//
// ======================================Cerebellum Bounds
template <typename PointStorageType, typename PointSetType>
bool
InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIPLCerebellarPointTypes(const std::string lmrkfilename)
{
  //  std::cout << "Reading IPL Cerebellar Bounds as  Landmark File " <<
  // lmrkfilename << std::endl;

  // Delete all old landmarks
  // DEBUG:  NOTE we can not delete after being called from Talairach.bnd!!!
  //        this->erase( this->begin(  ), this->end(  ) );
  // Open the files
  FILE * tempfile = fopen(lmrkfilename.c_str(), "r");

  if (tempfile == NULL)
  {
    return false;
  }
  constexpr unsigned short int FILE_BUFFER_SIZE = 256;
  char                         buffer[FILE_BUFFER_SIZE]; // Dummy Varible
  char *                       status = 0;
  // Read Header Information
  status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
  if (status == NULL)
  {
    return false;
  }
  if (strncmp(buffer, "IPL_HEADER_BEGIN", 16) != 0)
  {
    // Only files with a header of IPL_HEADER_BEGIN will be accepted
    return false;
  }
  // Read dimensions of image that landmarks are associated with
  // skip a head to landmark definitions
  status = buffer;
  while ((strncmp(buffer, "TALAIRACH_PARAMETER_HEADER_BEGIN", 32) != 0) && (status != NULL))
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_BEGIN found" << std::endl;
    }
  }

  bool xisset = false;
  bool yisset = false;
  bool zisset = false;
  // process looking for dims
  status = buffer;
  while ((strncmp(buffer, "TALAIRACH_PARAMETER_HEADER_END", 30) != 0) && (status != NULL))
  {
    constexpr unsigned short int TALAIRACH_PARAMETER_NAME_SIZE = 40;
    char                         Dummy[TALAIRACH_PARAMETER_NAME_SIZE];
    // Look for x
    if (strncmp(buffer, "TALAIRACH_PARAMETER_X_SIZE", 26) == 0)
    {
      sscanf(buffer, "%s %hu", Dummy, &(ImageDims[0]));
      xisset = true;
    }
    // Look for y
    if (strncmp(buffer, "TALAIRACH_PARAMETER_Y_SIZE", 26) == 0)
    {
      sscanf(buffer, "%s %hu", Dummy, &(ImageDims[1]));
      yisset = true;
    }
    // Look for z
    if (strncmp(buffer, "TALAIRACH_PARAMETER_Z_SIZE", 26) == 0)
    {
      sscanf(buffer, "%s %hu", Dummy, &(ImageDims[2]));
      zisset = true;
    }
    // if(!(xisset && yisset && zisset))
    // {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    // }
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
  }

  ImageDims[3] = 1;
  if (!(xisset && yisset && zisset))
  {
    std::cout << "ERROR must set all TALAIRACH_PARAMETER_[XYZ]_SIZE values" << std::endl;
  }

  // Skip ahead to data
  status = buffer;
  while ((strncmp(buffer, "IPL_HEADER_END", 14) != 0) && (status != NULL))
  {
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before IPL_HEADER_END found" << std::endl;
    }
  }

  // Actually read the data
  // First read the SLA point
  {
    PointType SLAPoint; // Superior Left Anterior
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("SLAPoint", buffer, ImageDims, ImageRes, SLAPoint) == false)
    {
      return false;
    }
    {                     // NOTE: for Cerebellum AC_Point is meaningless
      PointType AC_Point; // Antieror Commisure
      status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
      if (status == NULL)
      {
        std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
      }
      // if(process_bnd_point("AC_Point", buffer, ImageDims, ImageRes,
      // AC_Point)==false) {return false;}
    }

    PointType FV_Point; // Posterior Commisure
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("FV_Point", buffer, ImageDims, ImageRes, FV_Point) == false)
    {
      return false;
    }

    PointType IRPPoint; // Inferior Right Posterior
    status = fgets(buffer, FILE_BUFFER_SIZE, tempfile);
    if (status == NULL)
    {
      std::cout << "ERROR: End of file reached before TALAIRACH_PARAMETER_HEADER_END found" << std::endl;
    }
    if (process_bnd_point("IRPPoint", buffer, ImageDims, ImageRes, IRPPoint) == false)
    {
      return false;
    }
    // Process to extend 4 points to all possible combinations
    PointType TempPnt = SLAPoint;
    {
      std::string CurrentLandmarkName;
      //                                            R/L
      //                             S/I                           A/P
      // Next 4 may conflict with PC Point and center plane
      // CurrentLandmarkName="Cbl_SLAPoint";    TempPnt[0]=SLAPoint[0];
      // TempPnt[1]=SLAPoint[1]; TempPnt[2]=SLAPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      // CurrentLandmarkName="Cbl_SLPPoint";    TempPnt[0]=SLAPoint[0];
      // TempPnt[1]=SLAPoint[1]; TempPnt[2]=IRPPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      // CurrentLandmarkName="Cbl_SRAPoint";    TempPnt[0]=IRPPoint[0];
      // TempPnt[1]=SLAPoint[1]; TempPnt[2]=SLAPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      // CurrentLandmarkName="Cbl_SRPPoint";    TempPnt[0]=IRPPoint[0];
      // TempPnt[1]=SLAPoint[1]; TempPnt[2]=IRPPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      CurrentLandmarkName = "Cbl_S_FV_APoint";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_S_FV_PPoint";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_SL_FV_Point";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_SR_FV_Point";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_S_FV_Point";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = SLAPoint[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;

      CurrentLandmarkName = "Cbl_FV_LAPoint";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_LPPoint";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_RAPoint";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_RPPoint";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_FV_APoint";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_FV_PPoint";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_L_FV_Point";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_R_FV_Point";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_FV_Point";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = FV_Point[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;

      // Next 4 may conflict with inferior plane
      // CurrentLandmarkName="Cbl_ILAPoint";    TempPnt[0]=SLAPoint[0];
      // TempPnt[1]=IRPPoint[1]; TempPnt[2]=SLAPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      // CurrentLandmarkName="Cbl_ILPPoint";    TempPnt[0]=SLAPoint[0];
      // TempPnt[1]=IRPPoint[1]; TempPnt[2]=IRPPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      // CurrentLandmarkName="Cbl_IRAPoint";    TempPnt[0]=IRPPoint[0];
      // TempPnt[1]=IRPPoint[1]; TempPnt[2]=SLAPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      // CurrentLandmarkName="Cbl_IRPPoint";    TempPnt[0]=IRPPoint[0];
      // TempPnt[1]=IRPPoint[1]; TempPnt[2]=IRPPoint[2];
      // (*this)[CurrentLandmarkName]=TempPnt;
      CurrentLandmarkName = "Cbl_I_FV_APoint";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = SLAPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_I_FV_PPoint";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = IRPPoint[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_IL_FV_Point";
      TempPnt[0] = SLAPoint[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_IR_FV_Point";
      TempPnt[0] = IRPPoint[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
      CurrentLandmarkName = "Cbl_I_FV_Point";
      TempPnt[0] = FV_Point[0];
      TempPnt[1] = IRPPoint[1];
      TempPnt[2] = FV_Point[2];
      (*this)[CurrentLandmarkName] = TempPnt;
    }
  }
  fclose(tempfile);
  if (this->size() > 0)
  {
    return true;
  }
  else
  {
    return false;
  }
}

template typename PointStorageType,
  typename PointSetType > bool
                          InverseConsistentLandmarks<PointStorageType, PointSetType>::ReadIPLCerebellarPointTypes(
                            const std::string lmrkfilename,
                            const int         XDim,
                            const int         YDim,
                            const int         ZDim,
                            const int         TDim)
{
  const bool testread = this->InverseConsistentLandmarks::ReadIPLPointTypes(lmrkfilename);

  this->InverseConsistentLandmarks::rescale(XDim, YDim, ZDim, TDim);
  return testread;
}

template <typename PointStorageType, typename PointSetType>
void
InverseConsistentLandmarks<PointStorageType, PointSetType>::AddExtendedPointTypes3D_UnitCube(
  const InverseConsistentLandmarks & input)
{
  // First create the center slice replicated in +-x and +-y directions
  InverseConsistentLandmarks centersliceLM = input;

  (*this) = input; // Used to reset image Dims
  bool is3d = false;
  for (typename InverseConsistentLandmarks::const_iterator iter = input.begin(); iter != input.end(); ++iter)
  {
    const PointStorageType xpoint = iter->second[0];
    const PointStorageType ypoint = iter->second[1];
    const PointStorageType zpoint = iter->second[2];
    if (zpoint > 1e-5) // This assumes image sizes less than 10000
                       // slices
    {
      is3d = true;
    }
    const PointStorageType tpoint = iter->second.GetT();
    const PointStorageType wpoint = iter->second.GetWeighting();
    std::string            base(iter->first);

    std::string pl; // pointlable
    // Center
    centersliceLM[base] = iter->second;
    pl = "Left";
    pl += base;
    // Left
    centersliceLM[pl] = PointType((xpoint - 1.0F), (ypoint - 0.0F), zpoint, tpoint, wpoint);
    pl = "TopLeft";
    pl += base;
    // TopLeft
    centersliceLM[pl] = PointType(xpoint - 1.0F, ypoint + 1.0F, zpoint, tpoint, wpoint);
    pl = "Top";
    pl += base;
    // Top
    centersliceLM[pl] = PointType(xpoint - 0.0F, ypoint + 1.0F, zpoint, tpoint, wpoint);
    pl = "TopRight";
    pl += base;
    // TopRight
    centersliceLM[pl] = PointType(xpoint + 1.0F, ypoint + 1.0F, zpoint, tpoint, wpoint);
    pl = "Right";
    pl += base;
    // Right
    centersliceLM[pl] = PointType(xpoint + 1.0F, ypoint - 0.0F, zpoint, tpoint, wpoint);
    pl = "BottomRight";
    pl += base;
    // BottomRight
    centersliceLM[pl] = PointType(xpoint + 1.0F, ypoint - 1.0F, zpoint, tpoint, wpoint);
    pl = "Bottom";
    pl += base;
    // Bottom
    centersliceLM[pl] = PointType(xpoint - 0.0F, ypoint - 1.0F, zpoint, tpoint, wpoint);
    pl = "BottomLeft";
    pl += base;
    // BottomLeft
    centersliceLM[pl] = PointType(xpoint - 1.0F, ypoint - 1.0F, zpoint, tpoint, wpoint);
  }
  for (typename InverseConsistentLandmarks::const_iterator iter = centersliceLM.begin(); iter != centersliceLM.end();
       ++iter)
  {
    const PointStorageType xpoint = iter->second[0];
    const PointStorageType ypoint = iter->second[1];
    const PointStorageType zpoint = iter->second[2];
    const PointStorageType tpoint = iter->second.GetT();
    const PointStorageType wpoint = iter->second.GetWeighting();
    std::string            base(iter->first);

    std::string pl;               // pointlable
    (*this)[base] = iter->second; // Center

    if (is3d) // Only Do this if it is a 3D image
    {
      pl = "Previous";
      pl += base;
      (*this)[pl] = PointType((xpoint), (ypoint), (zpoint - 1.0F), tpoint, wpoint);

      pl = "Next";
      pl += base;
      (*this)[pl] = PointType((xpoint), (ypoint), (zpoint + 1.0F), tpoint, wpoint);
    }
  }
  return;
}

template <typename PointStorageType, typename PointSetType>
void
InverseConsistentLandmarks<PointStorageType, PointSetType>::AddExtendedPointTypes3D_OnN(
  const InverseConsistentLandmarks & input,
  const int                          nx,
  const int                          ny,
  const int                          nz)
{
  // First create the center slice replicated in +-x and +-y directions
  InverseConsistentLandmarks centersliceLM = input;

  (*this) = input; // Used to reset image Dims
  bool is3d = false;
  for (typename InverseConsistentLandmarks::const_iterator iter = input.begin(); iter != input.end(); ++iter)
  {
    const PointStorageType xpoint = iter->second[0];
    const PointStorageType ypoint = iter->second[1];
    const PointStorageType zpoint = iter->second[2];
    if (zpoint > 1e-5) // This assumes image sizes less than 10000
                       // slices
    {
      is3d = true;
    }
    const PointStorageType tpoint = iter->second.GetT();
    const PointStorageType wpoint = iter->second.GetWeighting();
    std::string            base(iter->first);

    std::string pl; // pointlable
    // Center
    centersliceLM[base] = iter->second;
    pl = "Left";
    pl += base;
    // Left
    centersliceLM[pl] = PointType((xpoint - nx), (ypoint - 0.0F), zpoint, tpoint, wpoint);
    pl = "TopLeft";
    pl += base;
    // TopLeft
    centersliceLM[pl] = PointType(xpoint - nx, ypoint + ny, zpoint, tpoint, wpoint);
    pl = "Top";
    pl += base;
    // Top
    centersliceLM[pl] = PointType(xpoint - 0.0F, ypoint + ny, zpoint, tpoint, wpoint);
    pl = "TopRight";
    pl += base;
    // TopRight
    centersliceLM[pl] = PointType(xpoint + nx, ypoint + ny, zpoint, tpoint, wpoint);
    pl = "Right";
    pl += base;
    // Right
    centersliceLM[pl] = PointType(xpoint + nx, ypoint - 0.0F, zpoint, tpoint, wpoint);
    pl = "BottomRight";
    pl += base;
    // BottomRight
    centersliceLM[pl] = PointType(xpoint + nx, ypoint - ny, zpoint, tpoint, wpoint);
    pl = "Bottom";
    pl += base;
    // Bottom
    centersliceLM[pl] = PointType(xpoint - 0.0F, ypoint - ny, zpoint, tpoint, wpoint);
    pl = "BottomLeft";
    pl += base;
    // BottomLeft
    centersliceLM[pl] = PointType(xpoint - nx, ypoint - ny, zpoint, tpoint, wpoint);
  }
  for (typename InverseConsistentLandmarks::const_iterator iter = centersliceLM.begin(); iter != centersliceLM.end();
       ++iter)
  {
    const PointStorageType xpoint = iter->second[0];
    const PointStorageType ypoint = iter->second[1];
    const PointStorageType zpoint = iter->second[2];
    const PointStorageType tpoint = iter->second.GetT();
    const PointStorageType wpoint = iter->second.GetWeighting();
    std::string            base(iter->first);

    std::string pl;               // pointlable
    (*this)[base] = iter->second; // Center

    if (is3d) // Only Do this if it is a 3D image
    {
      pl = "Previous";
      pl += base;
      (*this)[pl] = PointType((xpoint), (ypoint), (zpoint - nz), tpoint, wpoint);

      pl = "Next";
      pl += base;
      (*this)[pl] = PointType((xpoint), (ypoint), (zpoint + nz), tpoint, wpoint);
    }
  }
  return;
}
} // namespace itk
#endif
